In [182]:
from netCDF4 import Dataset, num2date
import glob
from mpl_toolkits.basemap import Basemap
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import calendar
from scipy.stats import describe
from scipy import interpolate
from tqdm import tqdm_notebook as tqdm
import oceansdb
In [279]:
def load_netcdf(path = "Argo_South_60"):
    files = glob.glob("data/" + path + "/**/*.nc", recursive=True)

    lats = []
    lons = []
    datetimes = []
    sst = []
    max_depth = []

    for f in tqdm(files):
        d = Dataset(f)
        lat = d.variables["LATITUDE"][:]
        mask = lat < -60
        lon = d.variables["LONGITUDE"][:]
        lats.extend(lat[mask])
        lons.extend(lon[mask])
        juld = d.variables["JULD"][:]
        units = d.variables["JULD"].getncattr('units')
        dates = num2date(juld, units, "standard")
        datetimes.extend(dates)
        try:
            sst.extend(d.variables["TEMP_ADJUSTED"][mask,0])
        except:
            sst.extend(np.full(len(mask), np.nan))
        max_depth.extend(np.nanmax(d.variables["PRES_ADJUSTED"], axis=1))

    lats = np.array(lats)
    lons = np.array(lons)
    datetimes = np.array(datetimes)
    sst = np.array(sst)
    np.save(path + "_lats", lats)
    np.save(path + "_lons", lons)
    np.save(path + "_dts", datetimes)
    np.save(path + "_sst", sst)
    np.save(path + "_depth", max_depth)
    return lats, lons, datetimes, sst, max_depth

def plot(lats, lons, z = [], title = "Argo profiles south of 60S", cbtitle = "Number of points in bin", vmax=None):
    # setup north polar stereographic basemap.
    # The longitude lon_0 is at 6-o'clock, and the
    # latitude circle boundinglat is tangent to the edge
    # of the map at lon_0. Default value of lat_ts
    # (latitude of true scale) is pole.
    fig = plt.figure(figsize=(15,15))
    m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
    m.drawcoastlines()
    m.fillcontinents(color='black',lake_color='aqua')
    # draw parallels and meridians.
    m.drawparallels(np.arange(-60, 0, 20))
    m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
    #m.drawmapboundary(fill_color='aqua')
    plt.title("{} {}".format(len(lats), title))

    x, y = m(lons, lats)
    if len(z) == 0:
        hh, locx, locy = np.histogram2d(x, y, bins=100)
        # Sort the points by density, so that the densest points are plotted last
        z = np.array([hh[np.argmax(a<=locx[1:]),np.argmax(b<=locy[1:])] for a,b in zip(x,y)])
        #print(describe(z))
        idx = z.argsort()
        x, y, z = x[idx], y[idx], z[idx]
    #m.imshow(heatmap, interpolation='bicubic', cmap="jet")
    m.scatter(x, y, c=z, s=5, alpha=1, cmap="jet", vmax=vmax)
    cb = plt.colorbar()
    cb.set_label(cbtitle, rotation=270)
    plt.show()

def plot_time(dts, title, label):
    if type(dts) == tuple:
        dtmin = min(dts[0].min(), dts[1].min())
        dtmax = max(dts[0].max(), dts[1].max())
    else:
        dtmin = dts.min()
        dtmax = dts.max()
    days = (dtmax - dtmin).days
    fig, ax = plt.subplots(1, 1, figsize=(15,15))
    print(days)
    ax.hist(dts, bins=days, label=label)
    ax.xaxis.set_major_locator(mdates.YearLocator())
    ax.set_title(title)
    ax.legend()
    plt.show()
In [229]:
argo_lats, argo_lons, argo_dts, argo_sst, argo_depth = load_netcdf()
plot(argo_lats, argo_lons, vmax=100)
plot(argo_lats, argo_lons, argo_sst[:len(argo_lats)], title = "Argo SST", cbtitle = "Temperature °C")
plot(argo_lats, argo_lons, argo_depth[:len(argo_lats)], title = "Argo Depth", cbtitle = "Depth (decibar pressure)")
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:30: UserWarning: Warning: converting a masked element to nan.
/usr/local/lib/python3.6/dist-packages/numpy/core/numeric.py:591: UserWarning: Warning: converting a masked element to nan.
  return array(a, dtype, copy=False, order=order, subok=True)
In [269]:
print(min(argo_dts), max(argo_dts))
plot_time(argo_dts, "Argo float timestamps", "argo")
2001-12-21 15:53:45 2019-07-04 19:53:14
6404
In [260]:
seal_lats, seal_lons, seal_dts, seal_sst, seal_depth = load_netcdf("seal")
plot(seal_lats, seal_lons, title = "Seal data south of 60S", vmax=100)
plot(seal_lats, seal_lons, seal_sst[:len(seal_lats)], title = "Seal SST", cbtitle = "Temperature °C")
plot(seal_lats, seal_lons, seal_depth[:len(seal_lats)], title = "Seal Depth", cbtitle = "Depth (decibar pressure)")
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:30: UserWarning: Warning: converting a masked element to nan.
In [261]:
print(min(seal_dts), max(seal_dts))
plot_time(seal_dts, "Seal data timestamps", "seal")
2004-01-27 11:48:59.999999 2018-01-14 23:00:00.000003
5101
In [262]:
all_lats = np.concatenate((seal_lats, argo_lats))
all_lons = np.concatenate((seal_lons, argo_lons))
plot(all_lats, all_lons, title="Argo + seal data south of 60S", vmax=100)
In [280]:
plot_time((seal_dts, argo_dts), "Argo + seal data timestamps", ["seal", "argo"])
6404
In [142]:
def load_netcdf_grid(path = "Argo_South_60", resolution = 1, target_var = "temp", with_depth = False, filter_month = None):
    if with_depth:
        grid = np.full(shape=[15 * resolution, 360 * resolution, 200], fill_value=np.nan)
    else:
        grid = np.full(shape=[15 * resolution, 360 * resolution], fill_value=np.nan)
    files = glob.glob("data/{}/**/*.nc".format(path), recursive=True)

    for f in tqdm(files[:]):
        d = Dataset(f)
        lat = d.variables["LATITUDE"][:]
        mask = (lat < -60) & (lat > -74.5)
        if filter_month:
            juld = d.variables["JULD"][:]
            units = d.variables["JULD"].getncattr('units')
            dates = num2date(juld, units, "standard")
            datemask = np.array([d.month == filter_month for d in dates])
            mask &= datemask
        if any(mask):
            lat = np.round((np.abs(lat[mask]) - 60) * resolution).astype(int)
            lon = d.variables["LONGITUDE"][mask]
            lon = np.round((lon + 180) * resolution).astype(int)
            lon[lon == 360 * resolution] = 0
            if with_depth:
                pres = np.round(d.variables["PRES_ADJUSTED"][mask] / 10).astype(int)
                temp = d.variables["TEMP_ADJUSTED"][mask]
                if target_var == "sal":
                    try:
                        sal = d.variables["PSAL_ADJUSTED"][mask]
                    except:
                        print("No salinity for {}".format(f))
                        continue
            else:
                pres = d.variables["PRES_ADJUSTED"][mask]
                temp = d.variables["TEMP_ADJUSTED"][mask, 0]
            for x in np.unique(lon):
                for y in np.unique(lat):
                    if with_depth:
                        ptmask = (lon == x) & (lat == y)
                        for z in np.unique(pres[ptmask]):
                            if z>= 200 or np.ma.is_masked(z):
                                continue
                            depthmask = pres[ptmask] == z
                            if target_var == "temp":
                                mean_for_pt = np.ma.mean(temp[ptmask][depthmask])
                            elif target_var == "sal":
                                mean_for_pt = np.ma.mean(sal[ptmask][depthmask])
                            if not np.isnan(mean_for_pt) and not np.ma.is_masked(mean_for_pt):
                                grid[y, x, z] = np.nanmean((grid[y, x, z], mean_for_pt))
                    else:
                        ptmask = (lon == x) & (lat == y)
                        if target_var == "n_point":
                            v = np.sum(ptmask)
                            if not np.ma.is_masked(v):
                                grid[y, x] = np.nansum((grid[y, x], v))
                        elif target_var == "temp":
                            temps_at_pt = temp[ptmask]
                            v = np.ma.mean(temps_at_pt)
                            if not np.ma.is_masked(v):
                                grid[y, x] = np.nanmean((grid[y, x], v))
                        elif target_var == "depth":
                            pres_at_pt = pres[ptmask]
                            if len(pres_at_pt):
                                v = np.ma.max(pres[ptmask])
                                if not np.ma.is_masked(v):
                                    grid[y, x] = np.nanmax((grid[y, x], v))
    return grid

def plot_grid(grid, resolution = 1, title="Gridded SST mean for argo data", cbtitle="Temperature °C", vmin=None, vmax=None):
    fig = plt.figure(figsize=(15,15))
    m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
    m.drawcoastlines()
    m.fillcontinents(color='black',lake_color='aqua')
    # draw parallels and meridians.
    m.drawparallels(np.arange(-60, 0, 20))
    m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
    #m.drawmapboundary(fill_color='aqua')
    plt.title("{} at {}° resolution".format(title, 1/resolution))

    x = np.arange(-180, 180.1, 1/resolution)
    y = np.arange(-60, -75.1, 1/-resolution)
    x, y = np.meshgrid(x, y)
    x, y = m(x, y)

    m.pcolormesh(x, y, grid, cmap="jet", vmin=vmin, vmax=vmax)
    cb = plt.colorbar()
    cb.set_label(cbtitle, rotation=270)
    plt.show()

def interp_nans(grid, method='linear'):
    x = np.arange(0, grid.shape[1])
    y = np.arange(0, grid.shape[0])
    mask = np.isnan(grid)
    xx, yy = np.meshgrid(x, y)
    x1 = xx[~mask]
    y1 = yy[~mask]
    newarr = grid[~mask]
    interp = interpolate.griddata((x1, y1), newarr.ravel(), (xx, yy), method=method)
    if method == "linear":
        interp = interp_nans(interp, "nearest")
    return interp
In [9]:
resolution = 2
argo_pt_grid = load_netcdf_grid(target_var="n_point", with_depth = False, resolution = resolution)
seal_pt_grid = load_netcdf_grid(path="seal", target_var="n_point", with_depth = False, resolution = resolution)
In [10]:
combined_pt_grid = np.nansum((argo_pt_grid, seal_pt_grid), axis=0)
combined_pt_grid[np.isnan(argo_pt_grid) & np.isnan(seal_pt_grid)] = np.nan
In [61]:
plot_grid(argo_pt_grid, title="Number of points in gridded argo dataset", cbtitle="Number of points", vmax=50, resolution = resolution)
plot_grid(seal_pt_grid, title="Number of points in gridded seal dataset", cbtitle="Number of points", vmax=50, resolution = resolution)
plot_grid(combined_pt_grid, title="Number of points in gridded argo+seal datasets", cbtitle="Number of points", vmax=100, resolution = resolution)
In [222]:
resolution = 4
In [223]:
argo_temp_grid = load_netcdf_grid(resolution = resolution)
seal_temp_grid = load_netcdf_grid("seal", resolution = resolution)
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:59: RuntimeWarning: Mean of empty slice
In [224]:
combined_temp_grid = np.nanmean((argo_temp_grid, seal_temp_grid), axis=0)
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:1: RuntimeWarning: Mean of empty slice
  """Entry point for launching an IPython kernel.
In [225]:
plot_grid(argo_temp_grid, resolution = resolution)
plot_grid(seal_temp_grid, title = "Gridded SST mean for seal data", resolution = resolution)
plot_grid(combined_temp_grid, title = "Gridded SST mean for argo+seal data", resolution = resolution)
In [226]:
interp = interp_nans(combined_temp_grid)
In [227]:
plot_grid(interp, resolution = resolution, title="Interpolated gridded SST mean for argo + seal data")
In [155]:
argo_max_depth = load_netcdf_grid(target_var = "depth", resolution = resolution)
seal_max_depth = load_netcdf_grid("seal", target_var = "depth", resolution = resolution)
1350
748
In [159]:
plot_grid(argo_max_depth, resolution = resolution, title="Argo gridded max depth", cbtitle="Depth (dbar)")
plot_grid(seal_max_depth, resolution = resolution, title="Seal gridded max depth", cbtitle="Depth (dbar)", vmax=2000)
In [237]:
for month in range(1, 13):
    argo_temp = load_netcdf_grid(filter_month = month, with_depth=True)
    seal_temp = load_netcdf_grid("seal", filter_month = month, with_depth=True)
    combined_temp = np.nanmean((argo_temp, seal_temp), axis=0)
    np.save("argo_temp_grid_{}".format(month), argo_temp)
    np.save("seal_temp_grid_{}".format(month), seal_temp)
    np.save("combined_temp_grid_{}".format(month), combined_temp)
    interp_temp = interp_nans(combined_temp[:,:,0])
    argomean = np.round(np.nanmean(argo_temp), 2)
    sealmean = np.round(np.nanmean(seal_temp), 2)
    combmean = np.round(np.nanmean(combined_temp), 2)
    interpmean = np.round(np.nanmean(interp_temp), 2)
    title = "Argo 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], argomean)
    plot_grid(argo_temp[:,:,0], title=title)
    title = "Seal 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], sealmean)
    plot_grid(seal_temp[:,:,0], title=title)
    title = "Argo+Seal 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], combmean)
    plot_grid(combined_temp[:,:,0], title=title)
    title = "Interpolated Argo+Seal 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], interpmean)
    plot_grid(interp_temp, title=title)
1350
748
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:4: RuntimeWarning: Mean of empty slice
  after removing the cwd from sys.path.
1350
748
1350
748
1350
748
1350
748
1350
748
1350
748
1350
748
1350
748
1350
748
1350
748
1350
748
In [9]:
resolution = 2
try:
    argo_temp_grid_withdepth = np.load("argo_temp_grid_withdepth.npy")
    seal_temp_grid_withdepth = np.load("seal_temp_grid_withdepth.npy")
    combined_temp_grid_withdepth = np.load("combined_temp_grid_withdepth.npy")
except:
    argo_temp_grid_withdepth = load_netcdf_grid(with_depth=True, resolution = resolution)
    seal_temp_grid_withdepth = load_netcdf_grid("seal", with_depth=True, resolution = resolution)
    combined_temp_grid_withdepth = np.nanmean((argo_temp_grid_withdepth, seal_temp_grid_withdepth), axis=0)
    np.save("argo_temp_grid_withdepth", argo_temp_grid_withdepth)
    np.save("seal_temp_grid_withdepth", seal_temp_grid_withdepth)
    np.save("combined_temp_grid_withdepth", combined_temp_grid_withdepth)
    
In [91]:
dbars = [30, 50, 100, 150, 199]
vmax = 5
vmin = -3
for db in dbars:
    dbs = "from {}-{} dbar".format(db*10, (db+1)*10)
    plot_grid(argo_temp_grid_withdepth[:,:,db], title="Gridded argo temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
    plot_grid(seal_temp_grid_withdepth[:,:,db], title="Gridded seal temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
    plot_grid(combined_temp_grid_withdepth[:,:,db], title="Gridded seal+argo temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
    interp = interp_nans(combined_temp_grid_withdepth[:,:,db])
    plot_grid(interp, title="Interpolated gridded seal+argo temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
In [5]:
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
#m.drawcoastlines()
#m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, 0, 20))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
m.etopo()
plt.title("Bathymetry south of -60S")
Out[5]:
Text(0.5,1,'Bathymetry south of -60S')
In [220]:
def plot_cross_slope(grid, title="Temperature at 179° W", lon=-179, cbtitle="Temperature °C", resolution = 1):
    fig, ax = plt.subplots(1, 1, figsize=(20,10))
    y = np.arange(-60, -75, -1/resolution)
    z = np.arange(0, 2000, 10)
    levels = np.arange(np.nanmin(grid), np.nanmax(grid), .01)
    im = ax.contourf(y, z, grid.T, cmap="jet", levels = levels, extend="both")
    ax.set_xlabel("° Latitude")
    ax.set_ylabel("Depth (dbar)")
    ax.set_title("{} at {} resolution".format(title, 1 / resolution))
    cb = fig.colorbar(im)
    cb.set_label(cbtitle, rotation=270)
    
    db = oceansdb.ETOPO()
    y = np.arange(-60, -75, -.01)
    h = -db['topography'].extract(lat=y, lon=lon)["height"]
    #ax.plot(y, h)
    ax.fill_between(y, 5000, h, color='black')
    ax.set_ylim(2000, 0)
    ax.set_xlim(-60, -74.5)
    
    plt.show()

filtered = combined_temp_grid_withdepth[:,2,:]
interp = interp_nans(filtered)
plot_cross_slope(filtered, resolution = resolution)
plot_cross_slope(interp, title="Interpolated Temperature at 179° W", resolution = resolution)
In [94]:
try:
    argo_sal_grid_withdepth = np.load("argo_sal_grid_withdepth.npy")
    seal_sal_grid_withdepth = np.load("seal_sal_grid_withdepth.npy")
    combined_sal_grid_withdepth = np.load("combined_sal_grid_withdepth.npy")
except:
    argo_sal_grid_withdepth = load_netcdf_grid(with_depth=True, target_var="sal", resolution = resolution)
    seal_sal_grid_withdepth = load_netcdf_grid("seal", with_depth=True, target_var="sal", resolution = resolution)
    combined_sal_grid_withdepth = np.nanmean((argo_sal_grid_withdepth, seal_sal_grid_withdepth), axis=0)
    np.save("argo_sal_grid_withdepth", argo_sal_grid_withdepth)
    np.save("seal_sal_grid_withdepth", seal_sal_grid_withdepth)
    np.save("combined_sal_grid_withdepth", combined_sal_grid_withdepth)
In [143]:
dbars = [30, 50, 100, 150, 199]
vmax = 35
vmin = 33.5
for db in dbars:
    dbs = "from {}-{} dbar".format(db*10, (db+1)*10)
    interp = interp_nans(combined_sal_grid_withdepth[:,:,db])
    plot_grid(interp, title="Interpolated gridded seal+argo salinity data " + dbs, cbtitle="Salinity", resolution=resolution, vmin=vmin, vmax=vmax)
In [153]:
def dens_smow(T):
    '''
    Function calculates density of standard mean ocean water (pure water) using EOS 1980.
    INPUT: T = temperature [degree C (ITS-90)]
    OUTPUT: dens = density [kg/m^3]
    '''
    a0 = 999.842594;a1 = 6.793952e-2;a2 = -9.095290e-3;a3 = 1.001685e-4;a4 = -1.120083e-6;a5 = 6.536332e-9
    T68 = T * 1.00024
    dens = (a0 + (a1 + (a2 + (a3 + (a4 + a5*T68)*T68)*T68)*T68)*T68)
    return dens

def dens0(S, T):
    '''
    Function calculates density of sea water at atmospheric pressure
    USAGE:  dens0 = dens0(S,T)
    DESCRIPTION:
        Density of Sea Water at atmospheric pressure using
        UNESCO 1983 (EOS 1980) polynomial.
    INPUT:  (all must have same dimensions)
        S = salinity    [psu      (PSS-78)]
        T = temperature [degree C (ITS-90)]
    OUTPUT:
        dens0 = density  [kg/m^3] of salt water with properties S,T,
        P=0 (0 db gauge pressure)
    '''
    assert S.shape == T.shape
    T68 = T * 1.00024
    # UNESCO 1983 eqn(13) p17.
    b0 = 8.24493e-1;b1 = -4.0899e-3;b2 = 7.6438e-5;b3 = -8.2467e-7;b4 = 5.3875e-9
    c0 = -5.72466e-3;c1 = +1.0227e-4;c2 = -1.6546e-6
    d0 = 4.8314e-4
    dens = (dens_smow(T) + (b0+ (b1+(b2+(b3+b4*T68)*T68)*T68)*T68)*S + (c0+(c1+c2*T68)*T68)*S*np.sqrt(S)+ d0*(S**2) )
    return dens
In [166]:
dens = dens0(combined_sal_grid_withdepth, combined_temp_grid_withdepth)
dbars = [30, 50, 100, 150, 199]
vmax = 1028
vmin = 1026.5

for db in dbars:
    dbs = "from {}-{} dbar".format(db*10, (db+1)*10)
    interp = interp_nans(dens[:,:,db])
    print(np.nanmin(interp), np.nanmax(interp))
    plot_grid(interp, title="Interpolated gridded seal+argo density data " + dbs, cbtitle="Density", resolution=resolution, vmin=vmin, vmax=vmax)
1027.0181758808599 1028.1316975356945
1027.1160229324942 1028.0460193495496
1027.0661376067346 1027.902182726604
1027.1902021401745 1027.8602134606635
1027.0167057697147 1027.8479702251195
In [219]:
plot_cross_slope(dens[:,2,:], title="Density at 179° W", cbtitle="Density", resolution=resolution)
plot_cross_slope(interp_nans(dens[:,2,:]), title="Interpolated density at 179° W", cbtitle="Density", resolution=resolution)
In [ ]: